Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 February 2008 (MN M£X style file v2.2) 



Evolution of planetesimal discs and planets migration 



m 
O 
o 



> 

m 
O 
m 
o 

6 



13 



A. Del Popolo, 1,2,3 S. Ye§ilyurt 3 , N. Ercan 3 

1 Dipartimento di Matematica, Universita Statale di Bergamo, Piazza Rosate, 2-1 24129 Bergamo, ITALY 

2 Feza Giirsey Institute, P. O. Box 6 Qengelkoy, Istanbul, Turkey 

3 Bogazici University, Physics Department, 80815 Bebek, Istanbul, Turkey 

2 February 2008 



ABSTRACT 

In this paper, we further develop the model for the migration of planets introduced 
in Del Popolo et al. j^) and extended to time-dependent planetesimal accretion disks 
in Del Popolo & Ekcsi More precisely, the assumption of Del Popolo & Ekcsi 

fl(|) that the surface density in planetesimals is proportional to that of gas is released. 
Indeed, the evolution of the radial distribution of solids is governed by many processes: 
gas-solid coupling, coagulation, sedimentation, evaporation/condensation, so that the 
distribution of planetesimals emerging from a turbulent disk does not necessarily re- 



3). 



In order to describe this evolution 
4o|) which, using a series of sim- 



flect that of gas (e.g., Stepinski & Valageas 
we use a method developed in Stepinski & Valageas 
plifying assumptions, is able to simultaneously follow the evolution of gas and solid 
particles for up to 10 yr. This model is based on the premise that the transformation 
of solids from dust to planetesimals occurs through hierarchical coagulation. Then, 
the distribution of planetesimals obtained after 10 7 yr is used to study the migration 
rate of a giant planet through the migration model introduced in This allows 
us to investigate the dependence of the migration rate on the disk mass, on its time 
evolution and on the value of the dimcnsionlcss viscosity parameter a. We find that 
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in the case of disks having a total mass of 10~ 3 — 10 Mq, and 10~ 4 < a < 10 -1 , 
planets can migrate inward over a large distance while if < 10 -3 Mq the planets 
remain almost at their initial position for a > 10~ 3 and only in the case a < 10~ 3 
the planets move to a minimum value of orbital radius of ~ 2AU. Moreover, the ob- 
served distribution of planets in the period range 0-20 days can be easily obtained 
from our model. Therefore, dynamical friction between planets and the planetesimal 
disk provides a good mechanism to explain the properties of observed extra-solar giant 
planets. 

Key words: Planets and satellites: general; planetary system 



1 INTRODUCTION 



Marcy et al. Q 



The discovery of solar-like stars showing evidences for planets orbiting around them (Mayor & Queloz 12(1 . Marcy et 
Vogt et al. 4g, Butler et al. |6|) has greatly intensified the interest in understanding the formation and evolution of planetary 
systems, as well as the long-standing problem of the solar system origin. At the same time these discoveries have raised a 
number of questions about the formation mechanisms of such systems. Indeed, the extra-solar planets discovered so far are 
all more massive than Saturn, and most either orbit very close to their stars or travel on much more eccentric paths than any 
of the major planets in our Solar System. The present small sample can be broadly divided into three groups: 

a) Jupiter analogues, in terms of period P and semi-axis a, and with low eccentricities (such as 47 UMa); 

b) planets with highly eccentric orbits (such as 70 Vir); 

c) close-in giants (or 'hot Jupiters') within 0.1 AU whose orbits are largely circular (such as 51 Peg). 



The properties of these planets, most of which are Jupiter-mass objects, are difficult to explain using the standard model 



r 3 



for planet formation (Lissauer 24; Boss 4j). Current theories (Mizuno 



U 



2g| Bodenheimer & Pollack |3() predict that giant planets 



were formed by gas accretion onto massive (~ 15Me) rocky cores which themselves were the result of the accumulation of 



a large number of icy planetesimals. The most favorable conditions for this process are found beyond the so-called "snow 







n 



line" (Hayashi |ljj; Sasselov & Lecar l3ftf) . As a consequence, this standard model predicts nearly circular planetary orbits and 



giant planets distances > 1 AU from the central star where the temperature in the protostellar nebula is low enough for icy 



materials to condense (Boss 



Q; BossQ 



iQ) 



but see also Wuchterll52t Wuchterl l53l) . Thus, in the case of close-in giants, it is very 
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unlikely that such planets were formed at their present locations. Then, the most natural explanation for this paradox and 
for planets on very short orbits is that these planets have formed further away in the protoplanetary nebula and they have 
migrated afterwards to the small orbital distances at which they are observed. Some authors have also proposed scenarios 

n o 

in which migration and formation were concurrent <l4ol . We refer the reader to Del Popolo et al. |9j (hereafter DPI) and Del 
Popolo & Ekcsi^] (hereafter DP2) for a detailed discussion of the different mechanisms that have been proposed to explain 
the presence of planets at small orbital distances: (a) dynamical instabilities in a system of giant planets (Rasio & Ford 

(c) tidal interaction with a gaseous nebula (Goldreich & 



imicai 1 

'■'>■>: WVidcii.-ihilliijg <'v Mwznri •"">.]). (1>) ~i> iteration instability^ i 3(l . (c) ti< 
Tremaine Goldreich & Tremaine Q WardQ Lin et al.Q WardQ) 
a planetesimal disk (DPI) 



and (d) dynamical friction between the planet and 



In particular, in DPI and DP2 we showed that dynamical friction between a planet and a planetesimals disk is an 
important mechanism for planet migration and we pointed out that some advantages of the model are: 

a) Planet halt is naturally provided by the model. 

b) It can explain planets found at heliocentric distances of > 0.03 — 0.04 AU, or planets having larger values of eccentricity. 

c) It can explain metallicity enhancements observed in stars having planets in short-period orbits. 

d) Radial migration is possible with modest masses of planetesimal disks, in contrast with other models. 

For the planetesimal disk used in DPI, following Opik (33), we assumed that the surface density in planetesimals E s varies 
as E s (r) = Eq (lAU/r) 3//2 , where E©, the surface density at 1 AU, was a free parameter. In DP2 the previous assumption 
was substituted by a more reliable model for the disk, and in particular we used a time-dependent accretion disk, since it is 
widely accepted that the solar system at early phases in its evolution is well described by this kind of structure. An important 
assumption of DP2 was that the surface density in planetesimals remains proportional to that of gas: E s (r, t) oc E(r, t). 
However, it is well-known that the distribution of planetesimals emerging from a turbulent disk does not necessarily reflect 



that of gas (e.g., Stepinski & Valageas 



a3 



Stepinski & Valageas |40j) . Indeed, in addition to gas-solid coupling, the evolution 
of the distribution of solids is also governed by coagulation, sedimentation and evaporation / condensation. In order to take 
into account these effects we use the method developed in Stepinski & Valageas (40) which is able to simultaneously follow 
the evolution of gas and solid particles for up to 10 7 yr. The main approximation used in this model is to associate one grain 
size to a given radius and time. Then, we use the radial distribution of planetesimals predicted by this model to estimate the 
planet migration, which is computed as in Del Popolo et al. JiJ). 

This paper is organized as follows. In Sect l2.1l we describe the disk model we use to obtain the radial distribution of the 
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planetesimal disk reached after 10 7 yr. Then, in Sect l2.2l we briefly review the migration model introduced in Del Popolo et 
al. il. Finally, we describe our results in Sect [31 



2 DISK MODEL AND PLANET MIGRATION 



2.1 Gas and solids distribution and evolution 



It is well-known that protostellar disks around young stellar objects that have properties similar to that expected for the 
solar nebula are common: between 25% to 75% of young stellar objects in the Orion nebula seem to have disks with mass 
1(T 3 M Q < Aid < 10 _1 M Q and size 40 ± 20 AU flj). Moreover, observations of circumsteilar disks surrounding T Tauri stars 
support the view of disks having a limited life-span and characterized by continuous changes during their life. These evidences 
have led to a large consensus about the nebular origin of the Solar System. Moreover, it clearly appears necessary to model 
both the spatial and temporal changes of the disk (which cannot be handled by the minimum-mass model nor by steady-state 
models). Besides, one also needs to describe the global evolution of the solid material which constitutes, together with the 
gas, the protoplanetary disk. Of particular interest is the spatial distribution of material making up a planetary system, as 
this is about the only information the present observations, based on the Doppler technique, can provide. The knowledge of 
this distribution and its time evolution is important to understand how planets form and in this paper it is a key issue since 
we wish to study the planet migration due to the interaction between planets and the local distribution of solid matter. 

As usual, the time evolution of the surface density of the gas E is given by the familiar equation (e.g., Stepinski & 
Valageas 4fli : 



5E _ 3d_ 

dt r dr 



1/2 9 / x/2 

r dr V VtE > 







(1) 



where ft is the turbulent viscosity. Since v% is not an explicit function of time, but instead depends only on the local disk 
quantities, it can be expressed as v t = ft(E, r) and Eq.Q can be solved subject to boundary conditions on the inner and 
outer edges of the disk. The opacity law needed to compute u t is obtained from Ruden & Pollack <I36T) . Then, Eq.Q is solved 
by means of an implicit scheme. Note that the evolution of the gas is computed independently from the evolution of particles 
(which only make ~ 1% of the gas mass). Next, from E(r, t) we can algebraically find all other gas disk variables. 

Next, as described in Stepinski & Valageas (4^) the evolution of the surface density of solid particles E s is given by: 



dE s _ 3d_ 
dt r dr 



i/2 d 



r - — (v s T, s r 



1/2, 



dr 



+ 



ld_ 

r dr 



2rE s («0) s 



(2) 



The first diffusive term is similar to Eq.Q, where the effective viscosity v s is given by: 
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ft 



Vs = <k with Sc = ^ + Y 1 + y2 • ( 3 ) 

Here we introduced the Schmidt number Sc which measures the coupling of the dust to the gas turbulence. We also used 
the relative velocity v between a particle and the gas, the turbulent velocity K, the Keplerian angular velocity and the 
so-called stopping time t s . The dimensionless quantity (fikt s ) measures the coupling of the solid particles to the gas. Thus, 
small particles (size s < 1 mm) have fik^s <§C 1 and they are strongly coupled to gas. Therefore, they exhibit the same radial 
velocity. On the other hand, large particles (s > 10 4 cm) with Q^t B S> 1 show a much smaller radial drift. Finally, particles in 
the intermediate regime (s ~ 10 cm) with Q^t s ~ 1 exhibit a large inward radial velocity. Therefore, the evolution of the dust 
radial distribution can be significantly different from the behaviour of the gas, depending on the particle size (see Stepinski 
& Valageas ^ for a detailed study) . 

The second advective term in Eq. arises from the lack of pressure support for the dust disk as compared with the gas 
disk. Thus, it is proportional to the azimuthal velocity difference iJ^, between the dust and the gas. The average (..) s refers to 
the vertical averaging over the disk height weighted by the solid density. We refer the reader to Stepinski & Valageas (|4fll for 
a more detailed presentation, see also Kornet et al. 



by tl 
I J). 



In addition to the radial motion described by Eq.101, the dust surface density also evolves through evaporation/condensation 
and coagulation. In this article, following Stepinski & Valageas l 4^) we assume that the size distribution of solid particles at 
a given orbital radius and time is narrowly peaked around a mean value s(r,t). This is supported by numerical simulations 
( 29) which show that although a broad size distribution is maintained, most of the mass is always concentrated in the largest 
particles so that one can define a meaningful typical size s(r,t). Then, within this approximation coagulation does not in- 
fluence the dust surface density E s since it conserves the total mass of solids. Thus, the coagulation of solid particles only 
appears through the evolution of the radial distribution s(r,t) of the typical size of the dust grains. We model this process 



as in Stepinski & Valageas 



We must note that we only consider collisional coagulation and we disregard gravitational 
interactions which would come into play at late times when large planetesimals have formed. 

On the other hand, we take into account the evaporation of solid particles which takes place at the radius r eva p where 
the temperature reaches T eV ap- We also include the condensation of the vapor below T ovap onto the solid grains. The velocity 
of the vapor is equal to gas velocity but this component evolves in a specific way because of its own diffusion process and 
condensation. In this article we are mainly interested in the distribution of solids at small radii hence we consider only one 
species of solid particles: high-temperature silicates with T evap = 1350 K and a bulk density pbuik = 3.3 g cm" . Thus, in our 
simplified model we follow the evolution of three distinct fluids: the gas, the vapour of silicates and the solid particles. 
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In this fashion, we obtain the radial distribution of the planetesimal swarm after 10 7 yr. This yields the surface density 
of solids E s (r, t) and the mid-plane solid density p s (r,t). We also obtain the size distribution s(r,t) reached by hierarchical 
coagulation. Of course, at these late times where planetesimals have typically reached a size of a few km or larger, gravitational 
interactions should play a dominant role with respect to coagulation. However, if these interactions do not significantly affect 
the radial distribution of solids (note that the radial velocity of such large particles due to the interaction with the gas is 
negligible) we can still use the outcome of the fluid approach described above to study the migration of giant planets, as 
detailed below. 



2.2 Migration model 

In order to study the migration of giant planets we use the model developed in DPI (see also DP2). Since this model has 
already been described in these two papers, we only recall here the main points. We consider a planet revolving around a 
star of mass M* = IMq. As described in the introduction, in this paper we suppose that the formation mechanism is core 
accretion and we are then not interested by the scenarios that gravitational instability models could introduce in our model. 
The equation of motion of the planet can be written as: 



r = F + R 



(4) 



where the term Fq represents the force per unit mass from the star, while R is the dissipative force (i.e. the dynamical 
friction term-see Melita & Woolfson l27ll . If we assume a disk-shaped matter distribution with constant velocity dispersions 
<7|| (parallel to the plane) and crj_ (perpendicular to the plane) and with a ratio simply taken to be 2:1 (i.e. a^—2a±), then 

Q 

according to Chandrasekhar ( 7 ) and Binney (|2 ) we may write the force components as: 



F \\ = k \\ v n 

= B \\ V M\ 
F ± = k ± vj 

= B±v 1± 

where 

Bn = 



2\/2nnG 2 log Amiro2 (mi + m<2 



(5) 



2v27rnG log AmiW2 (mi + 7712 



dq 

() (l+g)2(l- e 2+ g)l/2 



(6) 



x exp 



"ill 1 



2<7n 1 + q 2af.l-e 2 + q 



(7) 



B± = 
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dq 



(l + g)(l-e2 + 9 )3/2 
'III 1 «?, 



2af,l + q 2af,l-e 2 + q 



(8) 



x cxp 
and 

e=(l-aiArf)°- s . (9) 

Here n is the average spatial density of field particles, mi is the mass of the test particle, m.2 is the mass of a field one, and 
log A is the Coulomb logarithm. Then, the frictional drag on the test particles may be written as: 

F = — fc||t)i||ei| — k±v 1 ±e± (10) 



i Q Ida & MakinoQ 



where ey and ei are two unit vectors parallel and perpendicular to the disk plane. 

Since the damping of eccentricity and inclination is more rapid than radial migration (Ida 155 Ida & Makino 16; DPI), 
we only deal with radial migration and we assume that the planet has negligible inclination and eccentricity (i p ~ e p ~ 0) 
and that the initial distance to the star of the planet is 5.2 AU. 1 For the objects lying in the plane, the dynamical drag is 
directed in the direction opposite to the motion of the particle and is given by: 

F ~ -fc||«||e||. (11) 

In order to calculate the effect of dynamical friction on the orbital evolution of the planet, we suppose that a\\=2a± 
and that the dispersion velocities are constant. If the planetesimals attain dynamical equilibrium, their equilibrium velocity 
dispersion, cr m , would be comparable to the surface escape velocity of the dominant bodies (Safronov 1969) such that 

where 9 is the Safronov number, m* and r* are the mass and radius of the largest planetesimals, (note that the planetesimals 
velocity dispersion, cr m , now introduced, is the velocity dispersion to be used for calculating the a which is present in the 
dynamical friction force). If instead we consider a two-component system, consisting of one protoplanet and many equal- 
mass planetesimals the velocity dispersion of planetesimals in the neighborhood of the protoplanet depends on the mass of 
the protoplanet. When the mass of the planet, M, is < 10 25 g, the value of < > 1 / 2 (being e m the eccentricity of the 
planetesimals) is independent of M therefore: 

e m ~ 2O(2m/3M ) 1/3 (13) 

1 The planet was set at an initial distance of 5.2 AU for similitudes with choices done in previous papers (e.g., Murray et al. 1998, 
Trilling et al. 
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(Ida & Makino 1993) where m is the mass of the planetesimals. When the mass of the planet reaches values larger than 
10 25_ 10 26 gatl AU] < e 2 > i/2 ig proportional to m 1/3 : 

e m ~ 6(M/3M ) 1/3 (14) 

(Ida & Makino 1993). As a consequence also the dispersion velocity in the disc is characterized by two regimes being it 
connected to the eccentricity by the equation: 

Cm ^ {eL + i m ) 1/2 «c (15) 

where i m is the inclination of planetesimals and v c is the Keplerian circular velocity. Following Stern and Del Popolo et 
al. (^) we assume that (i^) = (e^)/4. In the simulation we assume that the planetesimals have all equal masses, m, and that 
m « M, M being the planet mass. This assumption does not affect the results, since dynamical friction does not depend on 
the individual masses of these particles but on their overall density. Note that for m M the frictional drag F\\ obtained in 
Eq.0 does not depend on the mass m of the planetesimals since the velocity dispersion <r m only depends on the mass M of 
the giant planet while the explicit dependence on m of Eq.JSJ only involves the product nm — p s . Therefore, we do not need 
to follow the evolution of the size distribution of planetesimals. We merely use the planetesimal density p B reached after 10 7 
yr, assuming that the height of the planetesimal disk does not evolve significantly. 

An important point to discuss here is the back reaction of the planet on the swarm. As previously reported, in the 
calculation of the scattering of planetesimals by a protoplanet, Ida & Makino (1993) showed, by means of N-body simulations, 
that eccentricities, e m , inclinations, i m and velocity dispersions, <r m , of planetesimals in the vicinity of the protoplanet are 
strongly influenced by the mass of the protoplanet. In the early stage, random velocities of small planetesimals remain low 
during the growth of the protoplanet, since they are regulated by gravitational scatterings between planetesimals. When the 
protoplanet becomes massive enough (10 25 — 10 26 g) to influence the velocity distribution of small planetesimals, the random 
velocities, and velocity dispersion of planetesimals are heated by the protoplanet and become larger than in the early stage. 
Furthermore, the protoplanet would scatter neighboring planetesimals and give rise to a gap in the spatial distribution of 
planetesimals (see Fig. 3 of Ida & Makino 1993, and Fig. 1 of Tanaka & Ida 1997). As noticed by Rafikov (2001) the gaps seen 
in N-body simulations are never clean because random motion of planetesimals is naturally included, and this permits some of 
them to be present in the gap. The process of clearing a gap in a planetesimal disc around a massive body is analogous to gap 
formation in gaseous discs (Takeuchi et al. 1996; Rafikov 2001). The gap also reduce the growth rate of the protoplanet. In 
our calculation, this effect was taken into account, as previously reported. The width of the heated region is roughly given by 
4[(4/3)(e 2 n + i 2 n )a 2 + 12/i^a 2 ] 1 / 2 (Ida & Makino 1993) where a is the semi-major axis and h M = ( in^) 173 is the Hill radius 
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of the protoplanet. The effect on the drift velocities can be easily predicted observing that the increase in velocity dispersion 
of planetesimals around the protoplanet decreases the dynamical friction force (see Eq. HOB and consequently increases the 
migration time-scale. In order to perform an order of magnitude estimation, we suppose that the planet moves on a stable 
circular orbit. Using for simplicity Chandrasekhar (1943) formula, the drift velocity can be expressed as: 

dr 4ir log AG 2 M p(r)r 

(Chandrasekhar 1943, Binney & Tremaine 1987, Palmer et al. 1993, Hernandez & Gilmore 1998). Recalling the mass depen- 
dence of the velocity dispersion in the two regimes, we find that the drift velocity behaves as: 



dr J M M<W^g 

oc { (17) 

constant M > 10 25 g 



dt 



The linear dependence on the planet's mass, of the drift velocity, corresponds to the type I drift in the density wave approach 
(Ward 1997), while the part of the plot independent on the planet's mass corresponds to type II drift. The transition between 
the two regimes entails a velocity drop of between one to three orders of magnitude (according to the value of a). What 
previously described is shown in Fig. 1, where we show the drift velocity, 4jr, as function of mass, M, of the protoplanet for 
M<j = 0.1 and ct = 10~ 4 . As shown, objects having masses < have velocity drift increasing as M, while after a threshold 
mass any further mass increases begins to slow down the drift. As the threshold is exceeded the motion fairly abruptly 
converts to a slower mode in which the drift velocity is independent of mass. As previously explained, this behavior is due to 
the transition from a stage in which the dispersion velocity is independent of M to a stage in which it increases with M 1 ^ 3 
(Ida & Makino 1993). This last stage is known as the protoplanet-dominated stage. The phenomenon is equivalent to that 
predicted in the density wave approach (Goldreich & Tremaine 1980; Ward 1997). In this approach, the density wave torques 
repel material on either side of the protoplanet's orbit and attempt to open a gap in the disc. Only very large objects are 
able to open and sustain the gap. After gap formation, the drift rate of the planet is set by disc viscosity and is generally 
smaller than in absence of the gap. We stress that the decaying portion of the curve corresponding to the transition from the 
first to the second stage does not correspond to any particular model because following Ida & Makino (1993) we do not have 
information on the evolution of a in the transition regime. Even if not necessary, in order to have an independent check of 
the effects of back reaction, we performed another calculation. Suppose to have the following two situations: 

a) there is an initial gap in the disc; 

b) the planet is initially embedded in an unperturbed disc and must clear material in order to form a gap. 

According to Nelson et al. (1999) the quoted situations lead to slightly different results in migration, results that tend to 



10 A. Del Popolo, S. Ye§ilyurt & N. Ercan 




1000 



Figure 1. Drift velocity, ^j, as a function of mass. Velocities arc normalized to V = ^jj^ ^Mq (^) 3r ^ where is an Earth mass, £ 
the surface density, Q is the angular velocity and a the dispersion velocity. The assumed conditions are those considered appropriate for 
the Jovian region and assuming that Md = O.IMq, a = 10~ 4 . The line oc M correspond to type I drift described by Ward (1997) and 
its behavior is valid till ~ O.IM9 but beyond this value there is a transition to a behavior oc M° (type II motion) 



converge with increasing time (see Fig. 3 of Nelson et al. 1999). The small initial difference in migration rate is due to the fact 
that in situation b the clearing of the disc material leads to a period of more rapid migration. Using this result, we compared 
the migration obtained assuming: 

1) an initial gap as that shown in Fig. 9 of Nelson et al. (1999) and with a = constant (not given by Ida & Makino (1993) 
relation) with 

2) the migration obtained taking account of the velocity dispersion dependence above described. 
The results were in good agreement. 
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Figure 2. Evolution of the particle size in a disk of total gas mass jWj = O.OIMq with a = 0.01, panel (a), and a = 0.0001, panel (b), 
at t = 10 4 yr (dashed line), t = 10 6 yr (dotted line) and t = 10 7 yr (solid line). 



3 RESULTS 

3.1 General framework 

In this article, similarly to DPI and DP2, we are mainly interested in studying the planet migration due to the interaction with 
planetesimals. For this reason we assume that the gas is almost dissipated when the planet starts its migration. 2 Moreover 
we know that the behavior of gas and dust /planetesimals is different. Usually, the decline of gas mass near stars is more rapid 
than the decline of the mass of orbiting solid matter l5^). While the gas tends to be dissipated, (several evidences show that 

3(l) . the coagulation process induces an 



the disk lifetimes range from 10 s yr to 10 7 yr, see Strom et al. 1 12t Ruden & Pollack 
increase of the density of solid particles with time (see Figs l516l and gives rise to objects of increasing dimensions, as shown 
in Fig|U 

We clearly see that after times of the order of 10 6 — 10 7 yr, the coagulation process gives rise to large particles (> 10 5 cm) 
which have a small radial velocity, whence a negligible radial motion. This leads to a freezing of the solid surface density to 
the value reached at times of the order of 10 7 yr. In other terms, once solids are in the form of planetesimals, the gas coupling 



2 Clearly the effect of the presence of gas should be that of accelerating the loss of angular momentum of the planet and to reduce the 
migration time. In our case, there is still gas after 10 7 yrs, but it is in quantity inferior to that of planetesimals, expecially in the case 
of lower mass discs, in which it can be even two order of magnitudes less than planetesimals. Moreover, as noticed by Kominami & Ida 
(2002) dynamical friction and gravitational gas drag are essentially the same dynamical process and then the effect of the gravitational 
gas drag is incorporated in that of dynamical friction. Other forces, as aerodynamical gas drag can be neglected when compared to 
gravitational gas drag (Kominami & Ida 2002), in our case 
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Table 1. Properties of the initial gas disk 
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becomes unimportant and the radial distribution of solids does not change any more. This is why we do not need to calculate 
its evolution for times longer than 10 7 years. 3 Then the disk is populated by residual planetesimals for a longer period. Here 
it is important to stress that planetesimal formation is not independent from initial conditions. In particular, the final solid 
surface density depends in an intricate fashion on the initial disk mass Afd and on the turbulent viscosity parameter a, see 

n n 

Stepinski & Valageas 14(1) and Kornet et al. II19T) . 

In order to investigate the dependence of the giant planet migration on the properties of the protoplanetary disk we 
integrated the model introduced in the previous sections for several values of the initial disk surface density (i.e. several disk 

n 

masses), and different values of a. More precisely, as in Stepinski & Valageas (4pj) we consider an initial gas surface density 
of the form: 

Eo(r) = Ei [l + (r/n) 2 ] ~ 3 ' 78 + Ea^/IAU)" 1 ' 5 . (18) 

The quantities Ei, n and E2 are free parameters which we vary in order to study different disk masses. The values we use are 
given in Tab.l, where Md is the gas disk mass (in units of Mq), J50 is the disk angular momentum (in units of 10 50 g cm 2 s _1 ), 
Ei and E2 are in g cm -2 and rj is in AU. The first term in Eg. 1181 ensures that there is some mass up to large distances from 
the star, while the second term corresponds to the central concentration of the mass and sets the location of the evaporation 
radius. Note however that in any case the evaporation radius for the high-temperature silicates we study here remains of order 
0.1 AU. As explained in Sect. 2.1 we consider only one species of solids: high-temperature silicates with T eV a P = 1350 K. We 
initialize the dust subdisk at time t = 10 4 yr (i.e. after the gas distribution has relaxed towards a quasi-stationary state) by 
setting the solid surface density E s as: E s = 6 x 10 _3 E in order to account for cosmic abundance. 



3.2 Evolution of the gas disk 

We show in Fig|^]the evolution of the midplane gas density p for the disk of initial mass Md = 0.1M©, with four values of 
a from 10 _1 down to 10 -4 . This range covers all values of a that are conceivable in the context of protoplanetary disks. Of 

3 Note however that the size distribution of planetesimals keeps changing due to their mutual gravitational interaction. 
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course, we can check that the gas disk is rapidly depleted when the viscosity is high. We show in Fig0]the results we obtain 
for the gas in the case of a small mass disk: M d = 1O" 4 M . Obviously we recover the same overall trends. Note that the main 
difference between our low-mass and high-mass scenarios shows in the outer radius of the gas disk. Indeed, with the initial 
conditions defined in Tab.l massive disks are more extended (see the characteristic radius ri). On the other hand, the gas 
density below 0.5 AU does not change by much so that the evaporation radius is always of order 0.1 AU. 

It is to note that in Figs. 3b,c,d, and Fig. 4c the solid line that corresponds to the longest time is, at places, above the 
other lines. Since the gas evolution is governed by a diffusion equation, one should expect that the surface density decreases 
with time at any point in the disk. In reality it is not always true that the surface density must decrease with time at any 
point in the disk, although its evolution is indeed given by a diffusion equation. In fact, in the inner parts of the disk there is 
also a net mass flux inward (so that the gas disk looses some mass onto the star). Then, depending on the initial conditions 
and a, if there is a lot of mass at large radii when this matter diffuses down to a smaller radius it can lead to a local increase 
of the surface density. This is what happens in Fig. 3b. Most of the mass initially lies at 20-200 AU, then, this mass slowly 
diffuses both outward and inward. Thus, one can check that in this range of radii the density indeed decreases with time. 
However, at very large radii the density increases somewhat (more precisely the outer radius of the disk grows). On the other 
hand, at small radii the density may also increase locally (but this is not always the case) as a large amount of matter may 
flow inward. Finally, in the quoted figures, we plot the density which makes the discussion a bit more difficult since it depends 
on both E and H (which decreases somewhat for smaller E). However this should not change the basics of the arguments. 

3.3 Evolution of the solid subdisk 

Next, we show in Fig|K|and Fig[(j]the evolution of the solid midplane density for the two cases Md = 1O _1 M0 and Md = 
1O _4 M0, and four values of a. We can see a converged radial distribution of solids emerge at late times of order 10 6 yr. Note 
that although the radial distribution of planetesimals depends on the value of a, the total mass of solids in a disk is roughly 
independent of a since it remains approximately equal to the initial mass of solids in the disk. This means that solids are 
reshuffled within the disk but they are not lost into the star. 4 The value of a determines the radial distribution of solids: 
particles in a disk with a larger value of a (more turbulent disk) have larger inward radial velocities and consequently are 
locked into planetesimals closer to the star than particles in a less turbulent disk. Thus, the smaller the value of a the broader 

4 In fact, solids initially located close to the evaporation radius are lost but they constitute a small percentage of the total solid material, 
which is predominantly located in the outer disk. 
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Figure 3. (a) Evolution of gas midplane density p, for a disk with M<j = O.IMq and a = 0.1 at t = 10 4 yr (dashed line), t = 10 6 yr 
(dotted line) and t = 10 7 yr (solid line), (b) Same as Fig. 2a but with a = 0.01. (c) Same as Fig. 2a but with a = 0.001. (d) Same as 
Fig.2a but with a = 0.0001 



the final distribution of solids. The evaporation radius is located at ~ 0.1 AU while the outer limit of converged p s is about 
10 AU for a = 10 _1 and moves to ~ 70 AU for a = 10~ 4 , in the case M d = O.IMq. In the case M d = 0.0001M Q , the outer 
radius moves only from ~ 10 to ~ 20 AU, since the initial extension of the disk is smaller. 

The coagulation process gives rise to solids of 10 6 - 10 7 cm. In disks characterized by smaller values of a, and thus a 
more extended distribution of solids, the range of sizes goes from 10 6 — 10 7 cm at the evaporation radius down to 10 3 — 10 4 
cm at the outer limit (see FigHJ . This is because the coagulation process is less efficient at larger radii where the solid density 
is smaller and the velocity dispersion of the dust decreases (along with the gas temperature which governs the turbulent 
velocity). At later times these solids will continue to increase their sizes, but will not change their radial position, as they are 
already large enough to have a negligible radial motion. Note that the converged radial distribution of solids does not vary 
monotonically. There are bulges of matter near the evaporation radius as well as at the outer limit of p s (or E s ). These bulges 
are present in all cases, but become narrower for smaller values of a. Their existence is a consequence of an intricate and 
nonlinear interdependence between advection and coagulation, modulated by changing gas conditions and the character of 
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Figure 4. (a) Evolution of gas midplanc density p, for a disk with = O.OOOIMq and a = 0.1 at t = 10 4 yr (dashed line), t = 10 e yr 
(dotted line) and t = 10 7 yr (solid line), (b) Same as Fig. 3a but with a = 0.01. (c) Same as Fig.3a but with a = 0.001. (d) Same as 
Fig.3a but with a = 0.0001 



turbulence. The formation of the inner bulge also signals the "freezing" of the total mass of the solids in the disk, inasmuch as 
no more solids are lost to the vapor zone in subsequent disk evolution. They will either be captured by the bulge or come to 
rest by themselves at larger radii. Notice that it is a radial squeezing due to particle dynamics, rather than vertical squeezing 
due to sedimentation, that is primarily responsible for establishing the bulge and keeping particles from falling into the vapor 
zone. The same mechanism is responsible for the abrupt drop in p s , which we associate with the outer limit of solid matter 
distribution, as well as the presence of the p B bulge at the location of this drop. 

Summarizing, the results show that the distribution of matter in the disk is sensitive to the assumed initial conditions. 
After times shorter than 3.2 x 10 6 yr the radial distribution of planetesimal mass density p B settles and can only change on a 
much longer timescale by processes not considered in our model (like mutual gravitational interactions between planetesimals 
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Figure 5. (a) Evolution of solids midplane density p a for a disk with = O.IMq and a = 0.1 at t = 10 4 yr (dashed line), t = 10 6 yr 
(dotted line) and t = 10 7 yr (solid line), (b) Same as Fig. 4a but with a = 0.01. (c) Same as Fig. 4a but with a = 0.001. (d) Same as 
Fig.4a but with a = 0.0001 



proposed first by Safronov l37t) . 5 We refer the reader to Stepinski & Valageas 
discussions of the behavior of protoplanetary disks. 



and Kornet et al. 
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for more detailed 



3.4 Migration of a giant planet 

Now, using the converged radial distribution of planetesimals derived in the previous section, we display in Fig|7]the evolution 
of the semi-major axis a(t) of a 1 Mj planet in such a disk for a = 0.1 and Ma =0.1 (solid line), 0.01 (dotted line), 0.001 
(short-dashed line), O.OOOlAf© (long-dashed line), respectively. We recall here that the planet is initially located at 5.2 AU 
with i p ~ e p ~ 0. 

We clearly see that for a fixed value of a a disk of larger mass leads to a more rapid migration of the planet. This behavior 
is quite natural since a more massive disk obviously yields a stronger frictional drag. Thus, in the cases — 0.1,0.01Mq 

5 Note that 10 6 yr is the time required for S s (or p a ) to converge everywhere. However, this convergence is not uniform and can be 
achieved on a time scale as short as 10 4 yr in the innermost disk where the dust density is highest. 
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Figure 6. (a) Evolution of solids midplane density p s for a disc with = O.OOOIMq and a = 0.1 at t = 10 4 yr (dashed line), t = 10 6 yr 
(dotted line) and t = 10 7 yr (solid line), (b) Same as Fig. 5a but with a = 0.01. (c) Same as Fig. 5a but with a = 0.001. (d) Same as 
Fig.5a but with a = 0.0001 



the planet migrates to ~ 0.08 AU in ~ 1.5 x 10 9 yr and to ~ 0.03 AU in ~ 2.5 x 10 9 yr, respectively. When the planet 
arrives at this distance the dynamical friction switches off and its migration stops. The stopping is simply due to the inner 
radius of the planetesimal disk. The latter is set by the evaporation radius. Indeed, solid bodies cannot condense at such 
small orbital radii r <0.1 AU because the temperature is too high. Of course, this evaporation radius r ovap depends on the 
properties of the solid grains we consider. For instance, for ice particles we have r cvap ~ 1 AU (e.g., Stepinski & Valageas 4(| 
Kornet et al. llST> . In this article we wish to understand the small orbital radii of observed planets, over the range 0.03 — 0.15 
AU. Therefore, we are interested in the inner regions of the disk where only high-temperature silicates with T ovap ~ 1350 K 
survive. This is why we selected this component in this study. Then, the main point of Fig[7] is that the properties of solids 
known to exist in protoplanetary systems, together with reasonable density profiles for the disk, lead in a natural fashion to a 
characteristic radius in the range 0.03 — 0.2 AU for the final semi-major axis of the giant planet. Note in particular that this 
process naturally explains why the migration stops at such radii. 



For less massive disks, Md = 10 ,10 M©, the migration is slower and the planet has not enough time to migrate 
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below ~ 4 AU. Similarly to DPI and DP2, note that we plotted our calculation till 4.5 x 10 9 yr which is the typical age of 
protoplanetary disks like ours (the Sun is ~ 4.5 x 10 9 yr old). Moreover, the planetesimal disk should have cleared out by this 
time (Ida et al. [rl DPI). It is interesting to note that the planet moves closer to the central star in the case Md = O.OIM© 
than in the case Md = 0.1-Mq. This is due to the fact that less massive disks usually have a smaller surface density which leads 
to a smaller temperature. This in turn implies a smaller evaporation radius so that the planet can move closer to the star. 6 
In order to study the effect of viscosity on migration, we performed three other calculations with a = 10~ 2 , 10~ 3 and 10 -4 , 
also plotted in Fig|7| We can note that the dependence of the final radius on a is rather weak and non-monotonic because it 
competes with the dependence on the surface density whose precise value is an intricate function of a. On the other hand, 
the migration time usually increases going from a = 10 -4 up to a = 0.1. Indeed, as seen from Fig|^|and FigEJa smaller a 
can lead to a larger mid-plane solid density. This is partly due to the fact that the dust disk height is smaller because the 
turbulence measured by a is weaker. 

The results displayed in Fig[7|show that for a < 0.01 a Jupiter-like planet can migrate to a very small distance from the 
parent star, 0.03AU < r < 0.1AU, provided the disk mass is sufficiently large M<j > 1O _3 M0. Only in the cases M d <W~ 4 M , 
or Md < 1O _3 M0 with a>0.1, the interaction with the planetesimal disk is too weak to yield a significant migration. Note 
however that a giant Jupiter-like planet already makes a mass of order 10 _3 Mq. Hence such objects should arise from 
protoplanetary disks with a mass of the order of or larger than 10 -3 M©. That is the mass of the initial protoplanetary disk 
should be at least of the same order since we do not expect all the matter to end up within this giant planet during its 
formation stage. 7 Then, our results show that a final radius of 0.03AU < r < 0.1 AU is a natural outcome, provided the 
planetesimal disk has not been cleared off too early on (e.g., by gravitational scattering). 

Summarizing, the present model predicts that, unless the disk mass is very small Md ^ 10 -4 Mq, planets tend to move 
close to the parent star and to pile up to distances of the order of 0.03 — 0.04 AU. However, with some degree of fine-tuning 
it is also possible to find a planet at intermediate distances between its formation site and such small radii. 

Before concluding this section, we have to add some hints of the way migration was calculated in the case the disc mass, 

6 Indeed, we have the energy balance: T 4 oc Ez^tf^l where T e is the effective temperature, S the gas surface density, i>t the turbulent 
viscosity and the Keplerian angular velocity. Besides, the mid-plane temperature T c obeys the relation T 4 ~ tT 4 where r is the 
opacity. On the other hand, the opacity r is given by t ~ kE where k is the Rosseland opacity, while the turbulent viscosity scales as 
ut ~ aTc/Q,^, hence we obtain: T 3 cx kbE 2 . 

7 However, if the giant planet forms through a different process, e.g., the instability of a distinct cloud, which has no relation with the 
disk itself, then the latter may have any mass. 
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Md, is smaller than that of the protoplanet. One can think that less massive discs are not able to move a more massive planet. 
This idea is not completely correct, as we shall see. As we know, when the planet mass is less than or comparable to the 
local disc mass, the planet behaves as a representative particle in the disc, as shown by Lin & Papaloizou (1986). When the 
protoplanet has a mass larger than that of the disc the satellite acts as a dam against the viscous evolution of the disc, and 
can lead to a substantial change in the disc structure in the vicinity of the planet. The coupled disc-planet evolution in this 
case has been studied by Syer & Clarke (1995) and by Ivanov et al. (1999), who showed that the inertia of the planet plays 
an important role in this case (see Nelson et al. 1999, Eq. 9). The effect of the interaction between the secondary and the disc 
leads to accumulation of the disc matter in the region behind the protoplanet (see Eq. 12 and 39 of Syer & Clarke 1995). In 
order to calculate the migration, in the quoted case, we used the surface density given in Ivanov et al. (1999) in the model of 
migration described in the previous sections. 



3.5 Distribution of orbital periods 

Finally, we show in Fig[S]the predictions of our model for the distribution of planets in the inner part of the disk. More 
precisely, we plot the fraction of planets in the orbital period range 0-20 days calculated from the model described in the 
previous section, assuming a uniform probability distribution in the plane (log(a), log(Afd)) in the range 10 -4 < a< 10 _1 and 
10~ 4 < M d < 1O _1 M . This is a rather arbitrary choice but our point is simply to check whether the observed distribution can 
be explained by reasonable values for a and Md or whether it requires some fine tuning. The right panel represents the same 



We can see that we actually 



distribution obtained with the data given in www.exoplanets.org, see also Kuchner & Lecar 
obtain a reasonable agreement with the data. Of course, there is some degeneracy so that different probability distributions 
in the plane (log(a), log(A/d)) could yield the same results. However, the main point of Fig[S]is to show that the observed 
distribution of orbital periods can be easily recovered from our model, without any fine-tuning and with reasonable values 
for a and Md. Moreover, as can be seen from Fig0 within our model the peak at 3-4 days for the orbital period comes from 
disks with Afd ^ 1O _3 M0 and a < 10 -2 . As noticed in the previous section this result actually fits nicely with the data since 
we can expect Jupiter-like planets to form in such protoplanetary disks with Md>lO _3 M0 (note also that the dependence 
on a is rather weak). 

Of the ~ 20 planets with periods less than 20 days, 7 have peri ods in the range 3-4 days, and it seems that this trend 
is not an artifact of observational selection (see Kuchner & Lecar \2(\l . As stressed by Kuchner & Lecar (|2oJ), phenomena like 
the interactions among two planets and a star, that can leave a planet trapped by stellar tides into a circular orbit of ~ 0.04 
AU (Rasio & Ford I35T) . or halting of planet migration by loss of mass to the star (Trilling et al. |47i) , are rare. The migration 
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Figure 7. (a) The evolution of the semi-major axis a(t) of a Jupiter-mass planet, M = lMj, in a planetesimal disk for a = 0.1 and 
several values of Mj, 0.1 (solid line), 0.01 (dotted line), 0.001 (short-dashed line) and 0.0001Mq (long-dashed line), (b) Same as Fig. 6a 
but with a = 0.01. (c) Same as Fig. 6a but with a = 0.001. (d) Same as Fig.6a but with a = 0.0001. 
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through resonant interaction with planetesimals in the disk (Murray et al. 30) provides a natural way of halting the planet but 
unfortunately to change the orbit of a Jupiter-mass planet requires approximatively a Jupiter mass of planetesimals and a too 
massive disk. The other possibility to explain the observed distribution of planets seen in Figgis connected with a gas disk 
truncated at a temperature of 1500 K by the onset of Magneto-Rotational instability (Kuchner & Lecar 2^). In other terms, 
disk temperature determines the orbital radii of the innermost surviving planets. In our model, the distribution of internal 
planets is naturally connected to the evaporation of planetesimals at very short radii. As stressed in DPI, DP2, this model 
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Figure 8. (a) Fraction of planets having orbital periods in the range 0-20 days, calculated according the model of this paper, (b) Fraction 
of planets having orbital periods in the range 0-20 days, calculated according to data (see www.exoplanets.org). 



has not the drawback of Murray et al. <I3CI). n amely that of requiring a too large disk mass for migration and at the same 
time it has the advantage of Murray et al. (30) of having an intrinsic natural mechanism that provides halting of migration. 



4 CONCLUSIONS 

In this paper, we further developed the model for the migration of planets introduced in DPI and extended to time-dependent 
planetesimal accretion disks in DP2. After releasing the assumption of DP2 that the surface density of planetesimals is 
proportional to that of gas, we used a simplified model developed by Stepinski & Valageas (1996, 1997), that is able to 
simultaneously follow the evolution of gas and high-temperature silicates for up to 10 7 yr. Then we coupled this disk model to 
the migration model introduced in DPI in order to obtain the migration rate of the planets in the planetesimal disk and to 
study how the migration rate depends on the disk mass, on its time evolution and on the dimensionless viscosity parameter 
a. 

We found that in the case of disks having a total mass of Md > 1O _3 M0 planets can migrate inward over a large distance, 
while if Md < 1O _3 M0 the planets remain almost at their initial position. On the other hand, for Md ~ 1O -3 M0 a significant 
migration requires cv< 10 -2 . 

If the migration is efficient the planet usually ends up at a small radius in the range 0.03 — 0.1 AU which is simply set 
by the evaporation radius of the gaseous disk which gave rise at earlier times to the radial distribution of the planetesimal 
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swarm. Thus, our model provides a natural explanation for the small observed radii of extra-solar giant planets. In particular, 
the halting of the inward migration of the planet is intrinsic to this process and it does not require a second mechanism. 

Finally, we noticed that the observed distribution of planet periods in the range 0-20 days can be easily obtained within 
this framework without fine-tuning. In particular, such small radii naturally occur for Jupiter-like planets if the initial disk 
has a mass of the same order or larger, which is quite likely. In order to inhibit this process (so that Jupiter-like planets like 
our own remain at larger distances > 5 AU) the planetesimal disk must be cleared off over a time-scale of the order of or 
smaller than 10 9 yr (depending on the properties of the disk) or the disk mass must be rather small (i.e. smaller than lMj) 
which could suggest an alternative formation scenario for such a giant planet (i.e. not related with the disk itself) . 
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